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The density and the elastic stiffness coefficients of fee sohd 
argon at high pressures from 1 GPa up to 80 GPa are com- 
puted by first-principles pseudopotential method with plane- 
wave basis set and the generalized gradient approximation 
(GGA). The result is in good agreement with the experimen- 
tal result recently obtained with the Brillouin spectroscopy 
by Shimizu et al. [Phys. Rev. Lett. 86, 4568 (2001)]. The 
Cauchy condition was found to be strongly violated as in the 
experimental result, indicating large contribution from non- 
central many-body force. The present result has made it clear 
that the standard density functional method with periodic 
boundary conditions can be successfully applied for calculat- 
ing elastic properties of rare gas solids at high pressures in 
contrast to those at low pressures where dispersion forces are 
important. 
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Rare gas atoms are among the simplest substances in 
physical and chemical natures because of their closed- 
shell electronic configuration. Many physical properties 
of rare gases at low pressures have been predicted by us- 
ing ab initio or empirical two-body potentials such as the 
Lennard-Johnes potentials, which consist of short-range 
repulsive potential and long-range attractive dispersion 
potential ||l|,D. At low temperatures, the rare gas atoms 
form van der Waals' crystals of fee structure except he- 
lium, which crystallizes to hep structure. The crystal 
structure and the binding energy of rare gas solids have 
been determined very accurately by experiments. It has 
been found that two-body potentials can describe these 
natures very well. Though two-body potentials underes- 
timate the binding energy by few percents, inclusion of 
three-body potentials as a small correction reduces the 
error within one percent 

Recently, the development of the Brillouin spec- 
troscopy in conjunction with diamond-anvil cells (DAC) 
1^ has opened a door for investigating elastic properties 
of rare gas solids at high pressures, which might be impor- 
tant for earth and planetary sciences. The pressure de- 
pendence of the elastic stiffness coefficients has been ex- 
perimentally determined up to 33 GPa by Grimsditch at 
al. 1^ and up to 70 GPa by Shimizu et al. 0. They used 
the envelope method to determine the elastic stiffness co- 
efficients: The Brillouin frequency shifts are measured 
without identifying the crystal orientation. Then the 
acoustic velocities calculated from the frequency shifts 
scatter because the acoustic velocity in crystals depends 



on the propagating direction. The maximum and min- 
imum values of acoustic velocities at each pressure are 
determined and the envelope curves of the acoustic ve- 
locities are drawn. The elastic stiffness constants are 
determined by comparing the envelope curves with the 
solutions of elastic equation. Shimizu et al. 0] have also 
determined the elastic stiffness coefficients of solid argon 
from 1.6 GPa up to 4 GPa by using in situ Brillouin spec- 
troscopy 1^ , in which the crystal orientation is identified 
and the more accurate results are obtainable. From the 
results of the in situ Brillouin spectroscopy and the en- 
velope method at high pressures, Shimizu et al. |Q found 
that solid argon becomes harder than iron |^ and that 
the deviation from the Cauchy relation becomes signif- 
icant. The latter implies that the contribution of non- 
central many-body force becomes more and more impor- 
tant at higher pressures, and it cannot be treated any 
more as a small correction to the two-body potentials. 

However, by now, calculations of elastic properties at 
high pressures have been limited to those with empiri- 
cal two-body potentials ||,^,|l^ as at low pressures. In 
this paper, therefore, we study the elastic stiffness co- 
efficients of fee solid argon at high pressures by using 
first principles calculations with periodic boundary con- 
ditions |pT]-jl3|| at zero temperature and under constant 
pressures, which can treat the effects of many-body po- 
tentials of crystals in a simple and direct way. The va- 
lence wave functions are expanded in a plane wave basis 
set truncated at a kinetic energy of 560 eV . The electron- 
ion interactions are described by the Vanderbilt-type ul- 
trasoft pseudopotentials |jl^. The effects of exchange- 
correlation interaction are treated within the generalized 
gradient approximation of Perdew et al. (GGA-PBE) 
p5[ . The model consists of a fee unit cell containing 
four argon atoms. The Brillouin zones are sampled with 
8x8x8 Monkhorst-Pack k-points [|l6| by using time- 
reversal symmetry only. During the structural optimiza- 
tion, the enthalpy H = E -\- PV is minimized by varying 
the length of the lattice vectors, while the angles between 
the lattice vectors and the atomic positions in the unit 
cell are fixed. In the geometrical optimization, the to- 
tal stress tensor is reduced to the order of 0.001 GPa by 
using the finite basis-set corrections [|l^ . 

The elastic stiffness tensor Cijki relates the stress tensor 
Oij and strain tensor e^i by the Hooke's law. 



Cyfcie/ci {i,j,k,l = x,y,z). 
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Since the stress and strain tensors are symmetric, the 
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most general elastic stiffness tensor has only 21 non-zero 
independent components. For cubic crystals, they re- 



duce to three components, cn 



Cl2 



and 



C44 = c 

coefficients can be determined by computing the stress 
generated by forcing a small strain to the optimized unit 
cell 18, 1^]. The lattice vectors of the strained unit 
cell are determined from the lattice vectors of the op- 
timized unit cell by the relation a^ = (/ + e)ai, where / 
is the unit matrix and e is a strain tensor. Two strain 
tensors, 



(2) 

V / 

/ e/2 0\ 
= I e/2 I , (3) 





are used to determine the three elastic stiffness coeffi- 
cients, cii, C12, and C44 from Eq.(|l|), namely, = cue, 
(j^y = cr^ = Ci2e and a^y = 0446 where and 
are the stress resulting from the strain and e^, re- 
spectively. The and are calculated with e = 0, 
0.01 and 0.02 at each pressure, and fitted to a parabolic 
function of e to remove non-linear contributions. The 
convergence of the elastic stiffness coefficients with re- 
spect to the cut-off energy and the number of k-points 
were estimated of the order of 1 GPa by increasing the 
cut-off energy to 1120 eV and Monkhorst-Pack k-points 
to 12 X 12 X 12, respectively. 
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FIG. 1. p — P equations of state for fee solid argon. The 
solid eurve with elosed cireles represents the present result. 
Open circles represent the experimental result obtained by 
Shimizu et al. |7| with the Brillouin spectroscopy. 
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FIG. 2. The pressure dependence of the elastic stiffness 
coefficients, cn, C12, and C44, of fee solid argon. The solid lines 
with closed symbols represent the present result. The dashed 
lines represent experimental result at 295 K obtained with 
envelope method by Shimizu et al. Q . Open symbols indicate 
the result obtained with m situ Brillouin spectroscopy Q. 
Cross symbols indicate the result of self-consistent phonon 
calculation based on pair potentials by Grimsditch et al. Q|. 

Figure 2 shows the pressure dependence of the elastic 
stiffness coefficients, cn, C12, and C44. The elastic stiff- 
ness coefficients increase linearly with increasing pres- 
sure. These elastic stiffness coefficients satisfy the gen- 
eralized elastic stability criteria for cubic crystals under 
hydrostatic pressure p0|-p2|. 



Cn -I- 2ci2 > 0, C44 > 0, 



and 



cii - C12 > 0. (4) 



The agreement between the present and the experimental 
results for cn and C12 does not seem excellent between 10 
GPa and 70 GPa. However, the agreement looks better 
when the acoustic velocities are plotted as a function of 
the pressure as shown in Fig. 3. 



Figure 1 shows the p — P equations of state. The 
agreement between the experiment and the GGA cal- 
culation indicates that the lattice constant of solid argon 
are mainly determined by the balance between the short- 
range repulsive force and the external pressure. The van 
der Waals' force, which is not taken into account explic- 
itly in the GGA calculation, is negligible in this pressure 
range. 
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FIG. 3. The pressure dependence of the squares of acous- 
tic velocities, wiA.max, 'yiA,min> VrA.^a^. and -yrA.min, of fee 
solid argon. The solid curves with closed symbols represent 
the present result. Open symbols represent the experimental 
data obtained by Shimizu et al. Q. The envelope method 
was used for above 4 GPa and in situ Brillouin spectroscopy 
was used below 4 GPa. The dashed curves are the envelope 
curves for the experimental data [Q]. 

Figure 3 shows the pressure dependence of the squares 
of the acoustic velocities, which are related to the elastic 
stiffness coefficients by 



LA. max 
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V. 



TA,niin 



(Cii +2ci2+4c44)/(3p) 
Cll/P 

C44//3 

(Cll - Ci2)/(2p). 



(5) 
(6) 
(7) 
(8) 



The agreement between the present result and the ex- 
perimental one becomes better than in Fig. 2, though 
theoretical VTA,min above 10 GPa seems slightly smaller 
than the experimental data. The difference between the 
present result and the experimental result in Fig. 2 may 
depend on how the experimental envelope curves are 
drawn in Fig. 3. 
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FIG. 4. The pressure dependence of the elastic anisotropy, 
A — 2c44/(cii — C12), for fee solid argon. Open squares rep- 
resent the theoretical result by Grimsditch et al. Other 
symbols and lines has meaning as in Fig3. 

Figure 4 shows the pressure dependence of the elastic 
anisotropy A = 2C44/ (cn — C12), which is the ratio of two 
shear moduli C44 and (cn — ci2)/2, and which becomes 
unity for isotropic elasticity. The anisotropy calculated 
between 1.6 and 4 GPa is approximately three, and agrees 
well with the experimental result. Above 4 GPa, the ex- 
perimental anisotropy slowly decreases toward two, while 
the present result is almost constant up to 80 GPa. 
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FIG. 5. The pressure dependence of the deviation from 
the Cauchy relation, 5 = c\2 — C44 — 2P, for fee solid argon. 
Symbols and lines have meanings as in Fig. 4. 

The deviation from the Cauchy relation 5 = C12 — C44 — 
2P is a measure of the contribution from the non-central 
many-body force since the Cauchy relation C12 — C44+2P 
should be satisfied when interatomic potentials are purely 
central . Figure 5 shows the pressure dependence of 
6. The deviation 6 for the experimental result [Q becomes 
larger as the pressure increases, which indicates that non- 
central many-body force becomes more and more impor- 
tant at high pressures. The 6 of the present result agrees 
well with the experimental result, indicating that first 
principles calculation with plane wave basis set and pseu- 
dopotentials can correctly describe the many-body force. 
The S for the theoretical result of Grimsditch et al. is 
almost zero for all pressures since their theory is based 
on pair potentials. 

In summary, we computed the elastic properties of fee 
solid argon at high pressure by using first principles cal- 
culation at zero temperature with plane wave basis set, 
pseudopotentials and generalized gradient approximation 
for exchange-correlation interaction. We have shown that 
the standard density functional method with periodic 
boundary conditions at zero temperature and under con- 
stant pressures can be successfully applied for calculating 
elastic properties of rare gas solids at high pressures in 
contrast to those at low pressures where dispersion forces 
are important. Though the effects of thermal and zero- 
point vibrations were not included in the present calcu- 
lation, thermal effects are expected to be small at high 
pressures. The long-range van der Waals' force is not 
also taken into account explicitly. However, the lattice 
constant and the elastic properties of solid argon at high 
pressures are mainly determined by the balance between 
the short-range repulsive force and the external pressure, 
and the van der Waals' force will be negligible. 
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